library(tidyverse)
library(readstata13)
library(foreign)
library(readxl)
library(writexl)
library(lubridate)
library(PanelMatch)

setwd("/Volumes/Backup Plus/Storage-OrganizedFiles/Article1-ReplicationMaterials_FullReplication/3.RobustnessChecks")

################################################################
# Main Model: Full Results
################################################################

rm(list = ls())

load("OMC-ResultsMainModel_Updated.RData")

################################################################
# Preparing the data - Protests
################################################################
x1 <- summary(pe.results5.p2_maha5b)
x2 <- summary(pe.results5.p2_ps.m5)
x3 <- summary(pe.results5.p2_ps.w5)
x4 <- summary(pe.results5.p2_CBPSm5)
x5 <- summary(pe.results5.p2_CBPSw5)

data1a <- as.data.frame(x1$summary)
data1a <- data1a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1a <- data1a %>%
  mutate(time = c(0, 1, 2, 3, 4))

data2a <- as.data.frame(x2$summary)
data2a <- data2a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data2a <- data2a %>%
  mutate(time = c(0, 1, 2, 3, 4))

data3a <- as.data.frame(x3$summary)
data3a <- data3a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data3a <- data3a %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4a <- as.data.frame(x4$summary)
data4a <- data4a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4a <- data4a %>%
  mutate(time = c(0, 1, 2, 3, 4))


data5a <- as.data.frame(x5$summary)
data5a <- data5a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5a <- data5a %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Plots - Protest
################################################################

pdf(file = "Figure1-A.pdf",   
    width = 9, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))

plot(x = data1a$time,
     y = data1a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1a$time, x1 = data1a$time,
         y0 = data1a$lower_b, y1 = data1a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data2a$time,
     y = data2a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data2a$time, x1 = data2a$time,
         y0 = data2a$lower_b, y1 = data2a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data3a$time,
     y = data3a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data3a$time, x1 = data3a$time,
         y0 = data3a$lower_b, y1 = data3a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4a$time,
     y = data4a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data4a$time, x1 = data4a$time,
         y0 = data4a$lower_b, y1 = data4a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5a$time,
     y = data5a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data5a$time, x1 = data5a$time,
         y0 = data5a$lower_b, y1 = data5a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Estimated Effect of Protest", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Preparing the data - Violent Self-Protection
################################################################

x1b <- summary(pe.results5.v2_maha5b)
x2b <- summary(pe.results5.v2_ps.m5)
x3b <- summary(pe.results5.v2_ps.w5)
x4b <- summary(pe.results5.v2_CBPSm5)
x5b <- summary(pe.results5.v2_CBPSw5)

data1b <- as.data.frame(x1b$summary)
data1b <- data1b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1b <- data1b %>%
  mutate(time = c(0, 1, 2, 3, 4))

data2b <- as.data.frame(x2b$summary)
data2b <- data2b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data2b <- data2b %>%
  mutate(time = c(0, 1, 2, 3, 4))

data3b <- as.data.frame(x3b$summary)
data3b <- data3b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data3b <- data3b %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4b <- as.data.frame(x4b$summary)
data4b <- data4b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4b <- data4b %>%
  mutate(time = c(0, 1, 2, 3, 4))


data5b <- as.data.frame(x5b$summary)
data5b <- data5b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5b <- data5b %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Plots - VSP
################################################################
pdf(file = "Figure2-A.pdf",   
    width = 9, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))

plot(x = data1b$time,
     y = data1b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1b$time, x1 = data1b$time,
         y0 = data1b$lower_b, y1 = data1b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data2b$time,
     y = data2b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data2b$time, x1 = data2b$time,
         y0 = data2b$lower_b, y1 = data2b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data3b$time,
     y = data3b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data3b$time, x1 = data3b$time,
         y0 = data3b$lower_b, y1 = data3b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4b$time,
     y = data4b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data4b$time, x1 = data4b$time,
         y0 = data4b$lower_b, y1 = data4b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5b$time,
     y = data5b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data5b$time, x1 = data5b$time,
         y0 = data5b$lower_b, y1 = data5b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Estimated Effect of VSP", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Preparing the data - Sanctuary
################################################################

x1c <- summary(pe.results5.s2_maha5b)
x2c <- summary(pe.results5.s2_ps.m5)
x3c <- summary(pe.results5.s2_ps.w5)
x4c <- summary(pe.results5.s2_CBPSm5)
x5c <- summary(pe.results5.s2_CBPSw5)

data1c <- as.data.frame(x1c$summary)
data1c <- data1c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1c <- data1c %>%
  mutate(time = c(0, 1, 2, 3, 4))

data2c <- as.data.frame(x2c$summary)
data2c <- data2c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data2c <- data2c %>%
  mutate(time = c(0, 1, 2, 3, 4))

data3c <- as.data.frame(x3c$summary)
data3c <- data3c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data3c <- data3c %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4c <- as.data.frame(x4c$summary)
data4c <- data4c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4c <- data4c %>%
  mutate(time = c(0, 1, 2, 3, 4))


data5c <- as.data.frame(x5c$summary)
data5c <- data5c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5c <- data5c %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Plots - Sanctuary
################################################################

pdf(file = "Figure3-A.pdf",   
    width = 9, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))

plot(x = data1c$time,
     y = data1c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 10, 6),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1c$time, x1 = data1c$time,
         y0 = data1c$lower_b, y1 = data1c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data2c$time,
     y = data2c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data2c$time, x1 = data2c$time,
         y0 = data2c$lower_b, y1 = data2c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data3c$time,
     y = data3c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data3c$time, x1 = data3c$time,
         y0 = data3c$lower_b, y1 = data3c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4c$time,
     y = data4c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data4c$time, x1 = data4c$time,
         y0 = data4c$lower_b, y1 = data4c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5c$time,
     y = data5c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data5c$time, x1 = data5c$time,
         y0 = data5c$lower_b, y1 = data5c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Estimated Effect of Sanctuary", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Robustness Checks 1: Main Models with up to 10 Matches
################################################################

rm(list = ls())

load("RCh1-OMC_ResultsMainModel_10M.RData")

################################################################
# Preparing the data - Protests
################################################################
x1 <- summary(pe.results5.p2_maha5b)
x2 <- summary(pe.results5.p2_ps.m5)
x3 <- summary(pe.results5.p2_ps.w5)
x4 <- summary(pe.results5.p2_CBPSm5)
x5 <- summary(pe.results5.p2_CBPSw5)

data1a <- as.data.frame(x1$summary)
data1a <- data1a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1a <- data1a %>%
  mutate(time = c(0, 1, 2, 3, 4))

data2a <- as.data.frame(x2$summary)
data2a <- data2a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data2a <- data2a %>%
  mutate(time = c(0, 1, 2, 3, 4))

data3a <- as.data.frame(x3$summary)
data3a <- data3a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data3a <- data3a %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4a <- as.data.frame(x4$summary)
data4a <- data4a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4a <- data4a %>%
  mutate(time = c(0, 1, 2, 3, 4))


data5a <- as.data.frame(x5$summary)
data5a <- data5a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5a <- data5a %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Plots - Protest
################################################################
  
pdf(file = "Figure1-B.pdf",   
      width = 9, 
      height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))

plot(x = data1a$time,
     y = data1a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1a$time, x1 = data1a$time,
         y0 = data1a$lower_b, y1 = data1a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data2a$time,
     y = data2a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data2a$time, x1 = data2a$time,
         y0 = data2a$lower_b, y1 = data2a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data3a$time,
     y = data3a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data3a$time, x1 = data3a$time,
         y0 = data3a$lower_b, y1 = data3a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4a$time,
     y = data4a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data4a$time, x1 = data4a$time,
         y0 = data4a$lower_b, y1 = data4a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5a$time,
     y = data5a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data5a$time, x1 = data5a$time,
         y0 = data5a$lower_b, y1 = data5a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Estimated Effect of Protest", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Preparing the data - Violent Self-Protection
################################################################

x1b <- summary(pe.results5.v2_maha5b)
x2b <- summary(pe.results5.v2_ps.m5)
x3b <- summary(pe.results5.v2_ps.w5)
x4b <- summary(pe.results5.v2_CBPSm5)
x5b <- summary(pe.results5.v2_CBPSw5)

data1b <- as.data.frame(x1b$summary)
data1b <- data1b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1b <- data1b %>%
  mutate(time = c(0, 1, 2, 3, 4))

data2b <- as.data.frame(x2b$summary)
data2b <- data2b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data2b <- data2b %>%
  mutate(time = c(0, 1, 2, 3, 4))

data3b <- as.data.frame(x3b$summary)
data3b <- data3b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data3b <- data3b %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4b <- as.data.frame(x4b$summary)
data4b <- data4b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4b <- data4b %>%
  mutate(time = c(0, 1, 2, 3, 4))


data5b <- as.data.frame(x5b$summary)
data5b <- data5b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5b <- data5b %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Plots - VSP
################################################################
pdf(file = "Figure2-B.pdf",   
    width = 9, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))

plot(x = data1b$time,
     y = data1b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1b$time, x1 = data1b$time,
         y0 = data1b$lower_b, y1 = data1b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data2b$time,
     y = data2b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data2b$time, x1 = data2b$time,
         y0 = data2b$lower_b, y1 = data2b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data3b$time,
     y = data3b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data3b$time, x1 = data3b$time,
         y0 = data3b$lower_b, y1 = data3b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4b$time,
     y = data4b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data4b$time, x1 = data4b$time,
         y0 = data4b$lower_b, y1 = data4b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5b$time,
     y = data5b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data5b$time, x1 = data5b$time,
         y0 = data5b$lower_b, y1 = data5b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Estimated Effect of VSP", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Preparing the data - Sanctuary
################################################################

x1c <- summary(pe.results5.s2_maha5b)
x2c <- summary(pe.results5.s2_ps.m5)
x3c <- summary(pe.results5.s2_ps.w5)
x4c <- summary(pe.results5.s2_CBPSm5)
x5c <- summary(pe.results5.s2_CBPSw5)

data1c <- as.data.frame(x1c$summary)
data1c <- data1c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1c <- data1c %>%
  mutate(time = c(0, 1, 2, 3, 4))

data2c <- as.data.frame(x2c$summary)
data2c <- data2c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data2c <- data2c %>%
  mutate(time = c(0, 1, 2, 3, 4))

data3c <- as.data.frame(x3c$summary)
data3c <- data3c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data3c <- data3c %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4c <- as.data.frame(x4c$summary)
data4c <- data4c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4c <- data4c %>%
  mutate(time = c(0, 1, 2, 3, 4))


data5c <- as.data.frame(x5c$summary)
data5c <- data5c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5c <- data5c %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Plots - Sanctuary
################################################################

pdf(file = "Figure3-B.pdf",   
    width = 9, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))

plot(x = data1c$time,
     y = data1c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 10, 6),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1c$time, x1 = data1c$time,
         y0 = data1c$lower_b, y1 = data1c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data2c$time,
     y = data2c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data2c$time, x1 = data2c$time,
         y0 = data2c$lower_b, y1 = data2c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data3c$time,
     y = data3c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data3c$time, x1 = data3c$time,
         y0 = data3c$lower_b, y1 = data3c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4c$time,
     y = data4c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data4c$time, x1 = data4c$time,
         y0 = data4c$lower_b, y1 = data4c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5c$time,
     y = data5c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data5c$time, x1 = data5c$time,
         y0 = data5c$lower_b, y1 = data5c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Estimated Effect of Sanctuary", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Robustness Checks 2: Main Models with L = 4
################################################################

rm(list = ls())

load("RCh2-OMC_ResultsMainModel_L4.RData")

################################################################
# Preparing the data - Protests
################################################################
x1 <- summary(pe.results5.p2_maha5b)
x2 <- summary(pe.results5.p2_ps.m5)
x3 <- summary(pe.results5.p2_ps.w5)
x4 <- summary(pe.results5.p2_CBPSm5)
x5 <- summary(pe.results5.p2_CBPSw5)

data1a <- as.data.frame(x1$summary)
data1a <- data1a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1a <- data1a %>%
  mutate(time = c(0, 1, 2, 3, 4))

data2a <- as.data.frame(x2$summary)
data2a <- data2a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data2a <- data2a %>%
  mutate(time = c(0, 1, 2, 3, 4))

data3a <- as.data.frame(x3$summary)
data3a <- data3a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data3a <- data3a %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4a <- as.data.frame(x4$summary)
data4a <- data4a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4a <- data4a %>%
  mutate(time = c(0, 1, 2, 3, 4))


data5a <- as.data.frame(x5$summary)
data5a <- data5a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5a <- data5a %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Plots - Protest
################################################################

pdf(file = "Figure1-C.pdf",   
    width = 9, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))

plot(x = data1a$time,
     y = data1a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1a$time, x1 = data1a$time,
         y0 = data1a$lower_b, y1 = data1a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data2a$time,
     y = data2a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data2a$time, x1 = data2a$time,
         y0 = data2a$lower_b, y1 = data2a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data3a$time,
     y = data3a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data3a$time, x1 = data3a$time,
         y0 = data3a$lower_b, y1 = data3a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4a$time,
     y = data4a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data4a$time, x1 = data4a$time,
         y0 = data4a$lower_b, y1 = data4a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5a$time,
     y = data5a$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data5a$time, x1 = data5a$time,
         y0 = data5a$lower_b, y1 = data5a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Estimated Effect of Protest", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Preparing the data - Violent Self-Protection
################################################################

x1b <- summary(pe.results5.v2_maha5b)
x2b <- summary(pe.results5.v2_ps.m5)
x3b <- summary(pe.results5.v2_ps.w5)
x4b <- summary(pe.results5.v2_CBPSm5)
x5b <- summary(pe.results5.v2_CBPSw5)

data1b <- as.data.frame(x1b$summary)
data1b <- data1b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1b <- data1b %>%
  mutate(time = c(0, 1, 2, 3, 4))

data2b <- as.data.frame(x2b$summary)
data2b <- data2b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data2b <- data2b %>%
  mutate(time = c(0, 1, 2, 3, 4))

data3b <- as.data.frame(x3b$summary)
data3b <- data3b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data3b <- data3b %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4b <- as.data.frame(x4b$summary)
data4b <- data4b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4b <- data4b %>%
  mutate(time = c(0, 1, 2, 3, 4))


data5b <- as.data.frame(x5b$summary)
data5b <- data5b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5b <- data5b %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Plots - VSP
################################################################
pdf(file = "Figure2-C.pdf",   
    width = 9, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))

plot(x = data1b$time,
     y = data1b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1b$time, x1 = data1b$time,
         y0 = data1b$lower_b, y1 = data1b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data2b$time,
     y = data2b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data2b$time, x1 = data2b$time,
         y0 = data2b$lower_b, y1 = data2b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data3b$time,
     y = data3b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data3b$time, x1 = data3b$time,
         y0 = data3b$lower_b, y1 = data3b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4b$time,
     y = data4b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data4b$time, x1 = data4b$time,
         y0 = data4b$lower_b, y1 = data4b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5b$time,
     y = data5b$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data5b$time, x1 = data5b$time,
         y0 = data5b$lower_b, y1 = data5b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Estimated Effect of VSP", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Preparing the data - Sanctuary
################################################################

x1c <- summary(pe.results5.s2_maha5b)
x2c <- summary(pe.results5.s2_ps.m5)
x3c <- summary(pe.results5.s2_ps.w5)
x4c <- summary(pe.results5.s2_CBPSm5)
x5c <- summary(pe.results5.s2_CBPSw5)

data1c <- as.data.frame(x1c$summary)
data1c <- data1c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1c <- data1c %>%
  mutate(time = c(0, 1, 2, 3, 4))

data2c <- as.data.frame(x2c$summary)
data2c <- data2c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data2c <- data2c %>%
  mutate(time = c(0, 1, 2, 3, 4))

data3c <- as.data.frame(x3c$summary)
data3c <- data3c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data3c <- data3c %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4c <- as.data.frame(x4c$summary)
data4c <- data4c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4c <- data4c %>%
  mutate(time = c(0, 1, 2, 3, 4))


data5c <- as.data.frame(x5c$summary)
data5c <- data5c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5c <- data5c %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Plots - Sanctuary
################################################################

pdf(file = "Figure3-C.pdf",   
    width = 9, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))

plot(x = data1c$time,
     y = data1c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 10, 6),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1c$time, x1 = data1c$time,
         y0 = data1c$lower_b, y1 = data1c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data2c$time,
     y = data2c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data2c$time, x1 = data2c$time,
         y0 = data2c$lower_b, y1 = data2c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data3c$time,
     y = data3c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data3c$time, x1 = data3c$time,
         y0 = data3c$lower_b, y1 = data3c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4c$time,
     y = data4c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data4c$time, x1 = data4c$time,
         y0 = data4c$lower_b, y1 = data4c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5c$time,
     y = data5c$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data5c$time, x1 = data5c$time,
         y0 = data5c$lower_b, y1 = data5c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Estimated Effect of Sanctuary", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Robustness Checks 3: Civilians Killed by FARC (Observatorio de Derechos Humanos de la Presidencia)
################################################################

rm(list = ls())

load("RCh3-ODDHP_AlternativeModel_FARC.RData")

################################################################
# Preparing the data - Protests
################################################################
x1 <- summary(pe.results5.p2_maha5b)
x2 <- summary(pe.results5.p2_ps.m5)
x3 <- summary(pe.results5.p2_ps.w5)
x4 <- summary(pe.results5.p2_CBPSm5)
x5 <- summary(pe.results5.p2_CBPSw5)

data1a <- as.data.frame(x1$summary)
data1a <- data1a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1a <- data1a %>%
  mutate(time = c(0, 1, 2, 3, 4))

data2a <- as.data.frame(x2$summary)
data2a <- data2a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data2a <- data2a %>%
  mutate(time = c(0, 1, 2, 3, 4))

data3a <- as.data.frame(x3$summary)
data3a <- data3a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data3a <- data3a %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4a <- as.data.frame(x4$summary)
data4a <- data4a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4a <- data4a %>%
  mutate(time = c(0, 1, 2, 3, 4))


data5a <- as.data.frame(x5$summary)
data5a <- data5a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5a <- data5a %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Plots - Protest
################################################################

pdf(file = "Figure1-D.pdf",   
    width = 9, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))

plot(x = data1a$time,
     y = data1a$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-60, 60, 6),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1a$time, x1 = data1a$time,
         y0 = data1a$lower_b, y1 = data1a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data2a$time,
     y = data2a$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-60, 60, 6),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data2a$time, x1 = data2a$time,
         y0 = data2a$lower_b, y1 = data2a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data3a$time,
     y = data3a$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-60, 60, 6),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data3a$time, x1 = data3a$time,
         y0 = data3a$lower_b, y1 = data3a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4a$time,
     y = data4a$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-60, 60, 6),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data4a$time, x1 = data4a$time,
         y0 = data4a$lower_b, y1 = data4a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5a$time,
     y = data5a$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-60, 60, 6),
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data5a$time, x1 = data5a$time,
         y0 = data5a$lower_b, y1 = data5a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Estimated Effect of Protest", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Preparing the data - Violent Self-Protection
################################################################

x1b <- summary(pe.results5.v2_maha5b)
x2b <- summary(pe.results5.v2_ps.m5)
x3b <- summary(pe.results5.v2_ps.w5)
x4b <- summary(pe.results5.v2_CBPSm5)
x5b <- summary(pe.results5.v2_CBPSw5)

data1b <- as.data.frame(x1b$summary)
data1b <- data1b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1b <- data1b %>%
  mutate(time = c(0, 1, 2, 3, 4))

data2b <- as.data.frame(x2b$summary)
data2b <- data2b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data2b <- data2b %>%
  mutate(time = c(0, 1, 2, 3, 4))

data3b <- as.data.frame(x3b$summary)
data3b <- data3b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data3b <- data3b %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4b <- as.data.frame(x4b$summary)
data4b <- data4b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4b <- data4b %>%
  mutate(time = c(0, 1, 2, 3, 4))


data5b <- as.data.frame(x5b$summary)
data5b <- data5b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5b <- data5b %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Plots - VSP
################################################################
pdf(file = "Figure2-D.pdf",   
    width = 9, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))

plot(x = data1b$time,
     y = data1b$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-60, 60, 6),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1b$time, x1 = data1b$time,
         y0 = data1b$lower_b, y1 = data1b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data2b$time,
     y = data2b$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data2b$time, x1 = data2b$time,
         y0 = data2b$lower_b, y1 = data2b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data3b$time,
     y = data3b$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data3b$time, x1 = data3b$time,
         y0 = data3b$lower_b, y1 = data3b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4b$time,
     y = data4b$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data4b$time, x1 = data4b$time,
         y0 = data4b$lower_b, y1 = data4b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5b$time,
     y = data5b$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data5b$time, x1 = data5b$time,
         y0 = data5b$lower_b, y1 = data5b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Estimated Effect of VSP", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Preparing the data - Sanctuary
################################################################

x1c <- summary(pe.results5.s2_maha5b)
x2c <- summary(pe.results5.s2_ps.m5)
x3c <- summary(pe.results5.s2_ps.w5)
x4c <- summary(pe.results5.s2_CBPSm5)
x5c <- summary(pe.results5.s2_CBPSw5)

data1c <- as.data.frame(x1c$summary)
data1c <- data1c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1c <- data1c %>%
  mutate(time = c(0, 1, 2, 3, 4))

data2c <- as.data.frame(x2c$summary)
data2c <- data2c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data2c <- data2c %>%
  mutate(time = c(0, 1, 2, 3, 4))

data3c <- as.data.frame(x3c$summary)
data3c <- data3c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data3c <- data3c %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4c <- as.data.frame(x4c$summary)
data4c <- data4c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4c <- data4c %>%
  mutate(time = c(0, 1, 2, 3, 4))


data5c <- as.data.frame(x5c$summary)
data5c <- data5c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5c <- data5c %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Plots - Sanctuary
################################################################

pdf(file = "Figure3-D.pdf",   
    width = 9, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4,1.5, 0))
par(mfrow = c(1, 5))

plot(x = data1c$time,
     y = data1c$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-60, 60, 6),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1c$time, x1 = data1c$time,
         y0 = data1c$lower_b, y1 = data1c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data2c$time,
     y = data2c$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data2c$time, x1 = data2c$time,
         y0 = data2c$lower_b, y1 = data2c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data3c$time,
     y = data3c$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data3c$time, x1 = data3c$time,
         y0 = data3c$lower_b, y1 = data3c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4c$time,
     y = data4c$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data4c$time, x1 = data4c$time,
         y0 = data4c$lower_b, y1 = data4c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5c$time,
     y = data5c$estimate, 
     main = "",
     ylim = c(-60, 60), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), xlab = "", ylab = "",
     pch = 16, yaxt = "n", cex.axis=1.25)
segments(x0 = data5c$time, x1 = data5c$time,
         y0 = data5c$lower_b, y1 = data5c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

mtext("Mahalanobis", line = -0.2, at = 0.11, outer = TRUE, cex = 1)
mtext("PS \n Matching", line = -1.8, at = 0.31, outer = TRUE, cex = 1)
mtext("PS \n Weighting", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("CBPS \n Matching", line = -1.8, at = 0.71, outer = TRUE, cex = 1)
mtext("CBPS \n Weighting", line = -1.8, at = 0.91, outer = TRUE, cex = 1)

mtext("Estimated Effect of Sanctuary", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years relative to the administration of treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# End of this part
################################################################
